#ifndef __OCCVALUE_LAM_CU__
#define __OCCVALUE_LAM_CU__


#include"Settings.h"
#include"cutil_inline.h"
#include"cutil_math.h"
#include "Index.h"
// copied from DNACount.c

/*constant memory*/
__constant__ unsigned int CountValue[ALPHABET_SIZE + 1];
__constant__ int c_shiftTable[MAX_QUERY_LENGTH];
__constant__ BWT c_bwt[1];
__constant__ BWT c_revbwt[1];

/*texture memory*/
unsigned int* d_texOcc = NULL;
texture<unsigned int, 2, cudaReadModeElementType> texOcc;
texture<unsigned int, 2, cudaReadModeElementType> texRevOcc;


/*
__device__
char getCharFromQueries(unsigned int * Queries, unsigned int queryId, unsigned int pos)
{
	unsigned int shift, rowIndex;
	char c=0;

	rowIndex = pos / CHAR_PER_WORD;
	shift = BIT_PER_CHAR * ( CHAR_PER_WORD -1  - (pos % CHAR_PER_WORD));
	c = (( Queries[queryId * QUERY_SIZE_IN_WORD + rowIndex])>>shift)&3;

	return c;
}
*/

__device__
char getCharFromQueries(const unsigned int * Queries, unsigned int queryId, unsigned int pos)
{
	unsigned int shift, rowIndex;
	char c=0;

	rowIndex = pos / CHAR_PER_WORD;
	shift = BIT_PER_CHAR * ( CHAR_PER_WORD -1  - (pos % CHAR_PER_WORD));
	c = (( Queries[queryId * QUERY_SIZE_IN_WORD + rowIndex])>>shift)&3;

	return c;
}


/*Queries are just for "one" query*/
//TODO: make a lookup table for the shift
__device__
char getCharFromQueries(const unsigned int * Queries, unsigned int pos)
{
	//unsigned int shift = BIT_PER_CHAR*(CHAR_PER_WORD - 1  - (pos&(CHAR_PER_WORD - 1)));
	//return (( Queries[pos>>LOG2_CHAR_PER_WORD])>>shift)&3;
	
	volatile unsigned int temp = Queries[pos>>LOG2_CHAR_PER_WORD];
	return (temp>>c_shiftTable[pos])&3;

	//return (( Queries[pos>>LOG2_CHAR_PER_WORD])>>c_shiftTable[pos])&3;
}


__device__ unsigned int DNAOccCount_2Layer(unsigned int *dna, unsigned int index,
                        char character, unsigned int backward) {
        unsigned int wordToCount, charToCount;
        unsigned int i, sum=0;
        unsigned long long d, b, a, mask;
        wordToCount = index>>(5);
        charToCount = index - wordToCount * 32;
        dna -= backward * 8;
        unsigned long long *dd = ((unsigned long long *)dna);
        i = backward ? 3 : 0;
        while(wordToCount > 0) {
                d = dd[i];
                b = (d >> 1);
                b ^= ((character & 0x2) ? 0 : 0xFFFFFFFFFFFFFFFF);
                a = d ^ ((character & 0x1) ? 0 : 0xFFFFFFFFFFFFFFFF);
                a &= b;
                a &= 0x5555555555555555;
                sum += __popcll(a);
                i = (backward ? i-1 : i+1);
                wordToCount--;
        }
        if (charToCount > 0) {
                d = dd[i];
                b = (d >> 1);
                b ^= ((character & 0x2) ? 0 : 0xFFFFFFFFFFFFFFFF);
                a = d ^ ((character & 0x1) ? 0 : 0xFFFFFFFFFFFFFFFF);
                a &= b;
                if (charToCount < 16) {
                        mask = (0x00000000FFFFFFFF << (32-charToCount*2)) & 0x00000000FFFFFFFF;
                } else {
                        mask = (0xFFFFFFFF00000000 << (64-charToCount*2)) | 0x00000000FFFFFFFF;
                }
                if (backward)
                        mask = __brevll(mask);
                a &= (mask & 0x5555555555555555);
                sum += __popcll(a);
        }
        return sum;
}

__device__ unsigned int GPUDNAOccCount(unsigned int *dna, unsigned int index, char character, unsigned int backward) {
	unsigned int wordToCount, charToCount;
	unsigned int i;
	unsigned int sum = 0;
	unsigned long long d, b, a, mask;
	
	
	//wordToCount = index / 32;
	wordToCount = index>>(5);
	charToCount = index - wordToCount * 32;
	dna -= backward * 4;
	
	ulonglong2 dd;
	//ulonglong2 is a 16Byte = 128bit CUDA vector
	//.x is the first 8Byte while .y is the last.
#if BYTE_4
	unsigned long long dd1, dd2;
	dd1 = ((unsigned long long)dna[1] << 32 | dna[0]);
	dd2 = ( (unsigned long long)dna[3] << 32 | dna[2]);
	
	dd = make_ulonglong2(dd1, dd2);
#elif BYTE_16
	dd = *((ulonglong2*)dna);
#endif
	i = 0;
	if (wordToCount > 0) { // i = 0
		d = backward ? dd.y : dd.x;
		b = (d >> 1);
		b ^= ((character & 0x2) ? 0 : 0xFFFFFFFFFFFFFFFF);
		a = d ^ ((character & 0x1) ? 0 : 0xFFFFFFFFFFFFFFFF);
		a &= b;
		a &= 0x5555555555555555;
		sum += __popcll(a);
		i = 1;
	}
	
	if (wordToCount == 2 || charToCount > 0) {
		if (backward) i = 3 - i;
		if (wordToCount == 2) charToCount = 32;
		d = (i & 1) ? dd.y : dd.x;
		b = (d >> 1);
		b ^= ((character & 0x2) ? 0 : 0xFFFFFFFFFFFFFFFF);
		a = d ^ ((character & 0x1) ? 0 : 0xFFFFFFFFFFFFFFFF);
		a &= b;
		if (charToCount < 16) {
			mask = (0x00000000FFFFFFFF << (32-charToCount*2)) & 0x00000000FFFFFFFF;
			//mask = (0x00000000FFFFFFFF << (32-(charToCount<<1))) & 0x00000000FFFFFFFF;
		} else {
			mask = (0xFFFFFFFF00000000 << (64-charToCount*2)) | 0x00000000FFFFFFFF;
			//mask = (0xFFFFFFFF00000000 << (64-(charToCount<<1))) | 0x00000000FFFFFFFF;
		}
		if (backward)
			mask = __brevll(mask);
		a &= (mask & 0x5555555555555555);
		sum += __popcll(a);
	}
	return sum;
}

__device__ void DNAAllOccCount_2Layer(unsigned int *dna, unsigned int index,
                                        unsigned int backward, unsigned int occCount[]) {
        unsigned int wordToCount, charToCount;
        unsigned int i;
        unsigned long long d, b, a, mask;

        wordToCount = index >> 5;
        charToCount = index - wordToCount * 32;
        dna -= backward * 8;
        unsigned long long  *dd = ((unsigned long long*)dna);
        i = backward ? 3 : 0;
        while(wordToCount > 0) {
                d = dd[i];
                for(int j=0; j < ALPHABET_SIZE; j++) {
                        b = (d >> 1);
                        b ^= ((j & 0x2) ? 0 : 0xFFFFFFFFFFFFFFFF);
                        a = d ^ ((j & 0x1) ? 0 : 0xFFFFFFFFFFFFFFFF);
                        a &= b;
                        a &= 0x5555555555555555;
                        if (backward)
                                occCount[j] -= __popcll(a);
                        else
                                occCount[j] += __popcll(a);
                }
                i = backward ? i-1 :i+1;
                wordToCount--;
        }
        
	if (charToCount > 0) {
                d = dd[i];
                for(int j=0; j<ALPHABET_SIZE; j++) {
                        b = (d >> 1);
                        b ^= ((j & 0x2) ? 0 : 0xFFFFFFFFFFFFFFFF);
                        a = d ^ ((j & 0x1) ? 0 : 0xFFFFFFFFFFFFFFFF);
                        a &= b;
                        if (charToCount < 16) {
                                mask = (0x00000000FFFFFFFF << (32-charToCount*2)) & 0x00000000FFFFFFFF;
                        } else {
                                mask = (0xFFFFFFFF00000000 << (64-charToCount*2)) | 0x00000000FFFFFFFF;
                        }
                        if (backward)
                                mask = __brevll(mask);
                        a &= (mask & 0x5555555555555555);
                        if (backward)
                                occCount[j] -= __popcll(a);
                        else
                                occCount[j] += __popcll(a);
                }
        }
}


__device__ void GPUDNAAllOccCount(unsigned int *dna, unsigned int index, 
					unsigned int backward, unsigned int occCount[]) {
	unsigned int wordToCount, charToCount;
	unsigned int i;
	unsigned long d, b, a, mask;
	
	//wordToCount = index / 32;
	wordToCount = index >> 5;
	charToCount = index - wordToCount * 32;
	dna -= backward * 4;
	ulonglong2 dd;
#if BYTE_4
        unsigned long long dd1, dd2;
        dd1 = ((unsigned long long)dna[1] << 32 | dna[0]);
        dd2 = ( (unsigned long long)dna[3] << 32 | dna[2]);

        dd = make_ulonglong2(dd1, dd2);
#elif BYTE_16
        dd = *((ulonglong2*)dna);
#endif  
	
	i = 0;
	if (wordToCount > 0) { // i = 0
		d = backward ? dd.y : dd.x;
#pragma unroll 4
		for (int j = 0; j < ALPHABET_SIZE; ++j) {
			b = (d >> 1);
			b ^= ((j & 0x2) ? 0 : 0xFFFFFFFFFFFFFFFF);
			a = d ^ ((j & 0x1) ? 0 : 0xFFFFFFFFFFFFFFFF);
			a &= b;
			a &= 0x5555555555555555;
			if (backward)
				occCount[j] -= __popcll(a);
			else
				occCount[j] += __popcll(a);
		}
		i = 1;
	}
	
	if (wordToCount == 2 || charToCount > 0) {
		if (backward) i = 3 - i;
		if (wordToCount == 2) charToCount = 32;
		d = (i & 1) ? dd.y : dd.x;
#pragma unroll 4
		for (int j = 0; j < ALPHABET_SIZE; ++j) {
			b = (d >> 1);
			b ^= ((j & 0x2) ? 0 : 0xFFFFFFFFFFFFFFFF);
			a = d ^ ((j & 0x1) ? 0 : 0xFFFFFFFFFFFFFFFF);
			a &= b;
			if (charToCount < 16) {
				mask = (0x00000000FFFFFFFF << (32-charToCount*2)) & 0x00000000FFFFFFFF;
			} else {
				mask = (0xFFFFFFFF00000000 << (64-charToCount*2)) | 0x00000000FFFFFFFF;
			}
			if (backward)
				mask = __brevll(mask);
			a &= (mask & 0x5555555555555555);
			if (backward)
				occCount[j] -= __popcll(a);
			else
				occCount[j] += __popcll(a);
		}
	}
}
__device__ unsigned int OccMinorValue(unsigned int *occMinor, unsigned int MinorIndex, char c) {
        unsigned int MinorExplicitIndex = MinorIndex / 2;
        if(MinorIndex % 2) {
                return (occMinor[MinorExplicitIndex * ALPHABET_SIZE + c] & 0xFFFF);
        }
        else {
                return (occMinor[MinorExplicitIndex * ALPHABET_SIZE + c] >> 16);
        }
}

__device__ unsigned int BWTOccValue_2Layer(unsigned int *bwt, unsigned int *occMajor, unsigned int *occMinor,
                                        unsigned int index, char c, unsigned int inverseSa0) {
        index -= (index > inverseSa0);
        unsigned int occMajorIndex = (index + (OCC_MINOR_INTERVAL / 2) - 1)/OCC_MAJOR_INTERVAL;
        unsigned int occMinorIndex = ( index + (OCC_MINOR_INTERVAL / 2) - 1) / OCC_MINOR_INTERVAL;
        unsigned int occIndex = occMinorIndex * OCC_MINOR_INTERVAL;
        unsigned int occValue = occMajor[occMajorIndex * ALPHABET_SIZE + c] + OccMinorValue(occMinor, occMinorIndex, c);

        if(occIndex != index) {
                unsigned int cnt = DNAOccCount_2Layer(bwt + (occIndex>>LOG2_CHAR_PER_WORD),
                                                        __usad(index, occIndex, 0), c, occIndex > index);
                return occIndex > index ? occValue - cnt : occValue + cnt;
        } else {
                return occValue;
        }
}
// copied from BWT.c
__device__ unsigned int GPUBWTOccValue(unsigned int *bwt, unsigned int *occ, 
			unsigned int index, char c, unsigned int inverseSa0) {

	index -= (index > inverseSa0);
	//unsigned int occExplicitIndex = (index + GPU_OCC_INTERVAL_DIV_2 - 1) / GPU_OCC_INTERVAL;
	unsigned int occExplicitIndex = (index + GPU_OCC_INTERVAL_DIV_2 - 1)>>LOG2_GPU_OCC_INTERVAL;
	unsigned int occIndex = occExplicitIndex * GPU_OCC_INTERVAL;
#if TEXTURE
	/*using texture memory*/
	int elementId = occExplicitIndex * ALPHABET_SIZE + c ;
	unsigned int occValue = tex2D(texOcc, elementId&65535, elementId>>16);
#else
	/*using global memory*/
	unsigned int occValue = occ[occExplicitIndex * ALPHABET_SIZE + c];
#endif	
	
	if (occIndex != index) {
		//unsigned int cnt = GPUDNAOccCount(bwt + occIndex / CHAR_PER_WORD, 
		//					__usad(index, occIndex, 0), c, occIndex > index);
		unsigned int cnt = GPUDNAOccCount(bwt + (occIndex>>LOG2_CHAR_PER_WORD), 
							__usad(index, occIndex, 0), c, occIndex > index);
		return occIndex > index ? occValue - cnt : occValue + cnt;
	} else {
		return occValue;
	}
	
}

__device__ void AllOccMinorValue(unsigned int * occMinor, unsigned int MinorIndex, unsigned int * value) {
        unsigned int MinorExplicitIndex = MinorIndex / 2;
        if(MinorIndex % 2) {
                for(int j=0; j< ALPHABET_SIZE; j++) {
                        value[j] = (occMinor[MinorExplicitIndex * ALPHABET_SIZE + j] & 0xFFFF);
                }
        }
        else {
                for(int j=0; j< ALPHABET_SIZE; j++) {
                        value[j] = (occMinor[MinorExplicitIndex * ALPHABET_SIZE + j] >> 16);
                }
        }
}
__device__ void BWTAllOccValue_2Layer(unsigned int *bwt, unsigned int *occMajor, unsigned int *occMinor,
                                unsigned int index, unsigned int inverseSa0, unsigned int occCount[]) {
	index -= (index > inverseSa0);
        unsigned int occMajorIndex = (index + (OCC_MINOR_INTERVAL / 2) - 1)/OCC_MAJOR_INTERVAL;
        unsigned int occMinorIndex = ( index + (OCC_MINOR_INTERVAL / 2) - 1) / OCC_MINOR_INTERVAL;
        unsigned int occIndex = occMinorIndex * OCC_MINOR_INTERVAL;
        int i;
        
        AllOccMinorValue(occMinor, occMinorIndex, occCount);
        for(i=0; i<ALPHABET_SIZE; i++) {
                occCount[i] += occMajor[occMajorIndex * ALPHABET_SIZE + i];
        }

        if(occIndex != index) {
                DNAAllOccCount_2Layer(bwt + (occIndex >> LOG2_CHAR_PER_WORD),
                                        __usad(index, occIndex, 0), occIndex > index, occCount);
        }
}

__device__ void GPUBWTAllOccValue(unsigned int *bwt, unsigned int *occ, 
				unsigned int index, unsigned int inverseSa0, unsigned int occCount[]) {
	index -= (index > inverseSa0);
	//unsigned int occExplicitIndex = (index + GPU_OCC_INTERVAL_DIV_2 - 1) / GPU_OCC_INTERVAL;
	unsigned int occExplicitIndex = (index + GPU_OCC_INTERVAL_DIV_2 - 1)>>LOG2_GPU_OCC_INTERVAL;
	unsigned int occIndex = occExplicitIndex * GPU_OCC_INTERVAL;
#if TEXTURE
	int elementId = occExplicitIndex * ALPHABET_SIZE;
	for(int i=0; i<ALPHABET_SIZE; i++){
		occCount[i] = tex2D(texRevOcc, elementId&65535, elementId>>16);
		elementId++;
	}
#else
	*((uint4*)occCount) = *((uint4*)(occ + occExplicitIndex * ALPHABET_SIZE));
#endif	
	if (occIndex != index) {
		//GPUDNAAllOccCount(bwt + occIndex / CHAR_PER_WORD, 
		//				__usad(index, occIndex, 0), occIndex > index, occCount);
		GPUDNAAllOccCount(bwt + (occIndex>>LOG2_CHAR_PER_WORD), 
						__usad(index, occIndex, 0), occIndex > index, occCount);
	}	
}

	
__device__
bool bwtCalc(BWT *bwt, char baseId, unsigned int * start, unsigned int * end)
{
	if ((*start) > bwt->textLength || (*end) > bwt->textLength)
		return false;

//	unsigned int occ = CountValue[baseId];
	unsigned int occ = 0;
/*		
	*start = occ +  GPUBWTOccValue(bwt->bwtCode,bwt->GPUOccValue, (*start), baseId,bwt->inverseSa0) + 1;
	*end = occ + GPUBWTOccValue(bwt->bwtCode,bwt->GPUOccValue,((*end)+1), baseId,bwt->inverseSa0);
	*start = occ +  GPUBWTOccValue((*c_bwt).bwtCode, (*c_bwt).GPUOccValue, (*start), 
						baseId, (*c_bwt).inverseSa0) + 1;
	*end = occ + GPUBWTOccValue((*c_bwt).bwtCode, (*c_bwt).GPUOccValue, ((*end)+1), 
						baseId, (*c_bwt).inverseSa0);
*/
#if   LAYER_INDEX_1     
        *start = occ +  GPUBWTOccValue(bwt->bwtCode,bwt->GPUOccValue, (*start), baseId,bwt->inverseSa0) + 1;
        *end = occ + GPUBWTOccValue(bwt->bwtCode,bwt->GPUOccValue,((*end)+1), baseId,bwt->inverseSa0);
#endif

#if LAYER_INDEX_2
        *start = occ + BWTOccValue_2Layer(bwt->bwtCode, bwt->occMajorValue, bwt->occValue, (*start), baseId,
                                                bwt->inverseSa0) + 1;
        *end = occ + BWTOccValue_2Layer(bwt->bwtCode,bwt->occMajorValue,bwt->occValue, ((*end)+1),
                                                baseId, bwt->inverseSa0);
#endif			
					
	if (*start <= *end)
		return true;
	return false;
}


__device__
void bwtCalcAll(BWT *revbwt, unsigned int *start, unsigned int *end, unsigned int *occCount_start, 
		unsigned int * occCount_end, unsigned int * occCount)
{
	int k;
#if LAYER_INDEX_1       
                GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, (*start), revbwt->inverseSa0, occCount_start) ;
                GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, (*end)+1, revbwt->inverseSa0, occCount_end) ;
#endif

#if LAYER_INDEX_2
                BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, (*start),
                                revbwt->inverseSa0, occCount_start) ;
                BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, (*end) + 1,
                                revbwt->inverseSa0, occCount_end) ;
#endif
                occCount[ALPHABET_SIZE-1] = 0;
                for (k=ALPHABET_SIZE-2; k>=0; k--)
                {
                        occCount[k] = occCount[k+1] + (occCount_end[k+1] - occCount_start[k+1]);
                }
}

	__device__
bool bwtCalcWithOccValue(unsigned int *countValue,char baseId, unsigned int startOccValue,
					unsigned int endOccValue, unsigned int * start, unsigned int * end)
{
//	(*start) = CountValue[baseId] + startOccValue + 1;
//	(*end) = CountValue[baseId] + endOccValue;
	(*start) = startOccValue + 1;
	(*end) = endOccValue;

	if ((*start) <= (*end))
		return true;
	return false;
}


__device__ void ExactMatch(BWT *bwt, unsigned int *queries, unsigned int*saInterval, unsigned int queryId)
{
	unsigned int start, end;
	int i;
	unsigned int mtableSize= bwt->lookupTable->tableSize;
	char c;
	unsigned long long packedPattern = 0;
	saInterval[0] = 1;
	saInterval[1] = 0;

	for (i = 0; i < mtableSize ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(queries, queryId, QUERY_LENGTH-mtableSize+i);
		packedPattern |= (c & 3);
	}
	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];

	for (i=QUERY_LENGTH-mtableSize-1 ; i>=0 && start<=end; i--){
		c = getCharFromQueries(queries, queryId, i);
		bwtCalc(bwt, c, &start, &end);
	}
	if(start <= end){
		saInterval[0] = start;
		saInterval[1] = end;
	}
}
__device__ void ExactMatchOccTest(BWT *bwt, unsigned int *queries, unsigned int*saInterval, unsigned int queryId)
{
        unsigned int start, end;
        int i;
        unsigned int mtableSize= bwt->lookupTable->tableSize;
        char c;
        unsigned long long packedPattern = 0;

        for (i = 0; i < mtableSize ; i++) {
                packedPattern<<=2;
                c = getCharFromQueries(queries, queryId, QUERY_LENGTH-mtableSize+i);
                packedPattern |= (c & 3);
        }
        start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
        end = bwt->lookupTable->table[packedPattern];

	int j;
        for (i=QUERY_LENGTH-mtableSize-1 ; i>=0 && start<=end; i--){
		saInterval[j*2] = start;
		saInterval[j*2 + 1] = end;
		j++;
                c = getCharFromQueries(queries, queryId, i);
                bwtCalc(bwt, c, &start, &end);
        }
}

__global__ void bwtExactMatch(unsigned int *Queries, unsigned int queryNum, BWT *bwt, unsigned int* saInterval)
{
	int iteration=0, queryId, i;
	__shared__ unsigned int queries[QUERY_SIZE_IN_WORD * THREAD_SIZE];
	unsigned int interval[2];


	for(iteration =0; iteration < ceil(queryNum / (double)(gridDim.x * blockDim.x)); iteration++) {
		queryId = iteration * gridDim.x * blockDim.x + blockIdx.x * blockDim.x + threadIdx.x;
		if(queryId >= queryNum)
			break;
		for(i=0; i<QUERY_SIZE_IN_WORD; i++)
			queries[threadIdx.x * QUERY_SIZE_IN_WORD + i] = Queries[queryId * QUERY_SIZE_IN_WORD + i];

		ExactMatch(bwt, queries, interval, threadIdx.x);
		saInterval[queryId * 2] = interval[0];
		saInterval[queryId * 2 + 1] = interval[1];
		
//		ExactMatchOccTest(bwt, queries, saInterval, threadIdx.x);
	}
}

#endif /*__OCCVALUE_LAM_CU__*/

